A basic analysis of the impact of sample size, effect size, and outcome standard deviation on statistical power in a simple RCT.
Author
Jeff Milliman
Published
July 22, 2026
Introduction
In this post, I am going to run through doing a power analysis via simulation for a basic experimental design - i.e. with two groups (treatment and control) of equal size - in R. To see how covariate adjustment impacts power, I’ll run the analysis both adjusted and unadjusted. I’ll also vary the sample size, the treatment effect, and the standard deviation of the outcome distribution to see how changing these parameters impacts power. Finally, I’ll plot the exaggeration ratio -i.e. the Type M Error (Gelman and Carlin 2014) - to see how the combination of filtering on statistical significance and low statistical power produces estimates that are 2 to 3 times larger than the true effect.
What is statistical power?
My preferred definition of statistical power is the following: “the probability of getting a statistically significant result, conditional on a set of beliefs about the model” (Blair, Coppock, and Humphreys 2023). This avoids all of the confusion about Type 1 vs Type 2 error and controlling error rates. While there are mathematical formulas to calculate statistical power (see Rainey (2025)), my preferred approach to power analysis is via simulation.
In practical terms, this means that you vary a bunch of model parameters, such as sample size, covariates, the treatment effect, size of the treatment group, etc.., generate thousands of simulations for each scenario, and calculate the number of times that you detect a statistically significant effect for your independent variable of interest. Researchers often set a target of 80% power, which means that if you generated a 1000 simulations of a data generating process with a given effect, you would detect a statistically significant effect in 800 of those scenarios at your chosen alpha level - usually p < .05.
In general, statistical power tends to increase as:
the sample size increases
the treatment effect increases
the standard deviation of the distribution of your outcome variable decreases
relevant covariates are adjusted for in the analysis
Why is low power bad?
One reason researchers give for why low power is bad is that researchers want to be sure that they can actually detect an effect of interest before they run a study. If a study has low power then researchers can’t be sure whether there truly is no effect if the result is not statistically significant. As increasing the sample size to increase statistical power tends to reach a point of diminishing returns - i.e. the point in which increasing the sample size has a negligible impact on statistical power - researchers also want to make sure that they can allocate their resources efficiently to reduce costs.
A second reason is that low powered studies, when filtered only for those that are statistically significant, tend to produce exaggerated effect sizes - i.e. effect sizes that are larger than the true effect size (Gelman and Carlin 2014). This is bad because the academic literature tends to select for studies based on statistical significance, which in turn leads to overly optimistic expectations about the impact of policy interventions (DellaVigna and Linos 2022; Simpson 2023).
For more on statistical power, J-PAL and the Evidence in Governance and Policy Network (egap) have good resources on statistical power for social scientists. Links to J-PAL and egap resources on statistical power.
Basic DGP for simple experimental design
In my hypothetical research design, I will model a simple randomized experiment conducted on students with equally sized treatment and control groups. The intervention will be designed to raise student test scores from 3 to 10 points on a 100 point test. I’ll also adjust for three covariates: age, high income (binary indicator), and gender (female, binary indicator), vary the sample size from 30 to 100 by 10, and include a final sample of 200 observations.
Outcome Distribution Standard Deviation = 5, 10 (higher noise = lower power)
Treatment Effect (tr) = 3, 5, 7, 10 (treatment group scores between 3 to 10 points higher on average)
Covariate distributions
Age: between 13 and 15 years (equal probability of age 13, 14, or 15)
Female: binary indicator (1 = female, 0 = male) with 60% probability female
High Income: binary indicator (1 = high income, 0 = low income) with 30% probability high income
Covariate Values
Age = 1 (a 1 year increase in age increases test scores by 1 point)
Female = 10 (Females score 10 points higher than Males on average)
High Income = 10 (high income students score 10 points higher than low income students on average)
Before running the simulation, I created a little sample DGP to see if the parameters produce reasonable looking test scores.
Code
#load in librarieslibrary(tidyverse)library(broom) #extracting pvalueslibrary(furrr) #running sims in parallel#sample DGP for power analysisset.seed(71)n <-50#outcome distributionx <-rnorm(n = n, mean =30, sd =5)#error terme <-rnorm(n = n, mean =0, sd =1)#treatment assignment - prob .5 = equal treatmentd <-rbinom(n =n, size =1, prob = .5)#age - between 13 and 15age <-sample(13:15, size = n, replace =TRUE)#gender - 60% femalefemale <-rbinom(n =n, size =1, prob = .6)#hi income - 30% hi incomehi_income <-rbinom(n =n, size =1, prob = .3)#treatment effecttr <-10#set the outcomey <- x + (age*1) + (female*10) + (hi_income*10) + (d*tr) + edata <-cbind(y, d, age, female, hi_income) |>as.data.frame()ols <-lm(y ~1+ d + age + female + hi_income, data = data)#histogram of the simulated test scoreshist(data$y)
Code
print(summary(ols))
Call:
lm(formula = y ~ 1 + d + age + female + hi_income, data = data)
Residuals:
Min 1Q Median 3Q Max
-11.312 -2.849 0.038 3.375 9.916
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.4329 12.3673 0.278 0.78261
d 10.0692 1.3946 7.220 4.81e-09 ***
age 2.7402 0.8737 3.136 0.00301 **
female 12.0786 1.3795 8.756 2.81e-11 ***
hi_income 10.9113 1.4833 7.356 3.03e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.7 on 45 degrees of freedom
Multiple R-squared: 0.804, Adjusted R-squared: 0.7866
F-statistic: 46.16 on 4 and 45 DF, p-value: 2.26e-15
The histogram of the test scores (y) looks ok. The test scores are fairly uniformly spread between 30 and 80, although the distribution is skewed slightly to the left, with a higher peak around test scores from 45 to 55. The test scores are maybe a little bit on the low side, but the distribution looks reasonable. With the exception of the intercept, the coefficients are all in the ball park of the covariate values that were simulated.
Function to run Power Analysis
To run the power analysis, I wrap the DGP in a function and return the coefficient values, along with a binary indicator of statistical significance for the treatment variable.
Code
#wrap DGP in a functionpower_function <-function(tr, n, prob_tr, mean, sd, age_eff, female_eff, income_eff, sig_level){#set outcome distribution x <-rnorm(n = n, mean = mean, sd = sd)#error term e <-rnorm(n = n, mean =0, sd =1)#set treatment assignment d <-rbinom(n = n, size =1, prob = prob_tr)#age - between 13 and 15 age <-sample(13:15, size = n, replace =TRUE)#gender - 60% female female <-rbinom(n =n, size =1, prob = .6)#hi income - 30% hi income hi_income <-rbinom(n =n, size =1, prob = .3)#set the outcome y <- x + (age*age_eff) + (female*female_eff) + (hi_income*income_eff) + (d*tr) + e#grab the data df <-cbind(y, d, age, female, hi_income) |>as.data.frame()#run model#unadjusted ols <-lm(y ~1+ d, data = df)#adjusted ols_adj <-lm(y ~1+ d + age + female + hi_income, data = df)#grab coefficents#unadjusted tr_unadjust <-coef(ols)[2]#adjusted intercept <-coef(ols_adj)[1] tr_coef <-coef(ols_adj)[2] age_coef <-coef(ols_adj)[3] female_coef <-coef(ols_adj)[4] income_coef <-coef(ols_adj)[5]#grab pvalue pvalue <- broom::tidy(ols)$p.value[2] pvalue_adjust <- broom::tidy(ols_adj)$p.value[2]#set true effect arguments true_treat_eff <- tr true_agg_eff <- age_eff true_female_eff <- female_eff true_income_eff <- income_eff#make into dataframe out <-tibble("intercept"= intercept,"control_mean"= mean,"sd"= sd,"treat_eff_unadj"= tr_unadjust,"treat_eff_adj"= tr_coef,"true_treat_eff"= true_treat_eff,"age_eff"= age_coef,"true_age_eff"= true_agg_eff,"female_eff"= female_coef,"true_female_eff"= true_female_eff,"income_eff"= income_coef,"true_income_eff"= true_income_eff, "sample_size"= n,"treat_pvalue"= pvalue,"treat_pvalue_adjust"= pvalue_adjust) |>#create statistical signifigance indicatorsmutate(stat_sig =case_when(treat_pvalue <= sig_level ~"yes",.default ="no"),stat_sig_adjust =case_when(treat_pvalue_adjust <= sig_level ~"yes",.default ="no"),sig_bin =case_when(stat_sig =="yes"~1,.default =0),sig_bin_adjust =case_when(stat_sig_adjust =="yes"~1,.default =0))return(out)}
Running the Simulations
Next, I vary my simulation parameters, creating a combination of all possible parameters with a 1000 simulations per scenario. To speed things up, I run the function with parallel processing from the furrr package. The simulation takes about 90 seconds to 2 minutes to run on my laptop.
Code
# Create the grid of all combinations of sim conditionssim_args <-expand.grid(tr =c(3,5,7,10),prob_tr = .5,age_eff =1,female_eff =10,income_eff =10,#sample size - i.e. nn =c(30,40,50,60,70,80,90,100,200),mean =30,#vary the standard deviationsd =c(5,10),sig_level = .05,#set 1000 simsnsims =1:1000) |>#drop nsims column for pmap function dplyr::select(-nsims)#wrap in quiet to capture warningsquiet_sim <- purrr::quietly(power_function)##set the plan for parallel processing - i.e. with 6 coresplan(multisession, workers =8)#start sim timetictoc::tic()#iterate through the sim arguments with the power sim functionsims <- furrr::future_pmap(sim_args, quiet_sim,#set the seed .options =furrr_options(seed =07192026),#set the progress bar.progress =TRUE)#end timetime_out <- tictoc::toc()
199.64 sec elapsed
Code
#calculate total sim timesim_time <- time_out$callback_msg#turn off the parallel processplan(sequential)##Extract our coefficients from the estimates into one data framesims_coefficients <- sims |>map("result") |> purrr::list_rbind(names_to ="sim_id") #grab our summary results - treatment effects and powersummary_results <- sims_coefficients |>#group by true effect size, sample, outcome sd dplyr::group_by(true_treat_eff, sample_size, sd) |> dplyr::summarise(obs_treat_eff =mean(treat_eff_unadj, na.rm =TRUE),obs_treat_eff_adj =mean(treat_eff_adj, na.rm =TRUE),power =mean(sig_bin, na.rm =TRUE),power_adj =mean(sig_bin_adjust, na.rm =TRUE)) |>ungroup()#calculate the m error/exaggeration ratiom_bias <- sims_coefficients |>filter(stat_sig =="yes") |>#group by model, sample size, and sim distribution dplyr::group_by(true_treat_eff, sample_size, sd) |> dplyr::summarise(obs_treat_eff =mean(treat_eff_unadj, na.rm =TRUE),true_treat_eff =mean(true_treat_eff,na.rm =TRUE),obs_treat_eff_adj =mean(treat_eff_adj, na.rm =TRUE),#get 95% confidence intervals around the effectseff_lo =quantile(treat_eff_unadj, .0275),eff_hi =quantile(treat_eff_unadj, .975),eff_lo_adj =quantile(treat_eff_adj, .0275),eff_hi_adj =quantile(treat_eff_adj, .975)) |>ungroup() |>#calculate m error (exaggeration ratio)mutate(mbias = obs_treat_eff/true_treat_eff,#calculate CI around mbiasmbias_lo = eff_lo/true_treat_eff,mbias_hi = eff_hi/true_treat_eff,mbias_adj = obs_treat_eff_adj/true_treat_eff,#calculate CI around mbias adjustedmbias_adj_lo = eff_lo_adj/true_treat_eff,mbias_adj_hi = eff_hi_adj/true_treat_eff)
Plot of Power Curve
Code
colors <-c("Adjusted"="blue", "Unadjusted"="red", "80% power threshold"="black")#power curvepower_curve <- summary_results |>#drop sample size 200filter(sample_size <=100) |>mutate(power = power*100,power_adj = power_adj*100,true_treat_eff =paste0("effect size ",true_treat_eff),true_treat_eff =factor(true_treat_eff, levels =c("effect size 3", "effect size 5","effect size 7","effect size 10")),sd =paste0("outcome sd ", sd),sd =factor(sd, levels =c("outcome sd 5", "outcome sd 10"))) |>ggplot(aes(x = sample_size)) +geom_line(aes(y = power, color ="Unadjusted")) +geom_line(aes(y = power_adj, color ="Adjusted")) +geom_hline(aes(yintercept =80, color ="80% power threshold"),lwd=1,linetype =2) +labs(title ="Statistical power by effect size and outcome standard deviation", y ="Power (percent statistically significant)",x ="Sample Size",color ="Legend") +scale_color_manual(values = colors,breaks =c("Adjusted", "Unadjusted", "80% power threshold")) +scale_y_continuous(breaks =seq(from =0, to =100, by =20), limits =c(0,110))+scale_x_continuous(breaks =seq(from =30, to =100, by =10), limits =c(15,110)) +theme_bw() +theme(plot.title =element_text(size =14),legend.title =element_text(size =14), legend.text =element_text(size =14),axis.title.x =element_text(size =14),axis.title.y =element_text(size =14),axis.text.x =element_text(size =12),axis.text.y =element_text(size =12),strip.text.x =element_text(size =14),plot.caption =element_text(size =12)) +facet_wrap(~ true_treat_eff + sd, ncol =2)power_curve
Above, I plot the relationship between sample size and statistical power for both the adjusted and unadjusted models for each combination of effect size and outcome standard deviation. As expected, as the sample size increases statistical power increases across all scenarios. Adjusting for covariates also increases statistical power, although the impact of covariate adjustment on statistical power is less noticeable in scenarios where the outcome variable has a higher standard deviation (10 sd vs 5 sd).
Finally, increasing the treatment effect also increases statistical power. At a simulated treatment effect of 3, 80% power is only barely reached with covariate adjustment, an outcome standard deviation of 5, and a sample size of 100. As the treatment effect increases, the adjusted and unadjusted models are more likely to reach 80% power and increasing the outcome standard deviation has less of an impact on power.
Plot of Type M Error (Exaggeration Ratio)
Code
#mbias plot#pivot the type m error, confidence intervals lo and confidence intervals hi and join togerthermbias_plot <- m_bias |>filter(sample_size <=100) |>select(true_treat_eff:sd, mbias,mbias_adj) |>mutate(sample_size =as.factor(sample_size)) |>pivot_longer(cols =c(mbias,mbias_adj),names_to ="bias_type",values_to ="bias") |>mutate("Model type"=case_when(str_detect(bias_type, "adj") ~"Adjusted",.default ="Unadjusted")) |>select(-bias_type)#pivot confidence intervals lombias_plot2 <- m_bias |>filter(sample_size <=100) |>select(true_treat_eff:sd, mbias_lo, mbias_adj_lo) |>mutate(sample_size =as.factor(sample_size)) |>pivot_longer(cols =c(mbias_lo, mbias_adj_lo),names_to ="bias_type",values_to ="mbias_lo") |>mutate("Model type"=case_when(str_detect(bias_type, "adj") ~"Adjusted",.default ="Unadjusted")) |>select(-bias_type)#pivot confidence intervals himbias_plot3 <- m_bias |>filter(sample_size <=100) |>select(true_treat_eff:sd, mbias_hi, mbias_adj_hi) |>mutate(sample_size =as.factor(sample_size)) |>pivot_longer(cols =c(mbias_hi, mbias_adj_hi),names_to ="bias_type",values_to ="mbias_hi") |>mutate("Model type"=case_when(str_detect(bias_type, "adj") ~"Adjusted",.default ="Unadjusted")) |>select(-bias_type)#creat join listjoin_list <-list(mbias_plot, mbias_plot2, mbias_plot3)#join into one dfmbias_plot_df <-reduce(join_list, left_join, by =c("true_treat_eff", "sample_size", "sd", "Model type"))#drop plot dfsrm(mbias_plot, mbias_plot2, mbias_plot3)#plot of type m error by effect size and outcome sdmbias_plot <- mbias_plot_df |>mutate(true_treat_eff =paste0("effect size ",true_treat_eff),true_treat_eff =factor(true_treat_eff, levels =c("effect size 3", "effect size 5","effect size 7","effect size 10")),sd =paste0("outcome sd ", sd),sd =factor(sd, levels =c("outcome sd 5", "outcome sd 10"))) |>ggplot(aes(x = sample_size, y = bias, group =`Model type`, colour =`Model type`)) +geom_point(size =2, position =position_dodge(width = .9)) +geom_errorbar(aes(ymin = mbias_lo, ymax = mbias_hi), position =position_dodge(width = .9), size = .7) +geom_hline(yintercept =1, linetype =2) +scale_y_continuous(breaks =seq(from =-3, to =6, by =1), limits =c(-3,6)) +geom_hline(aes(yintercept =1, color ="True effect = Observed effect"),lwd=1,linetype =2) +scale_color_manual(values =c("Adjusted"="blue", "Unadjusted"="red","True effect = Observed effect"="black"),breaks =c("Adjusted", "Unadjusted", "True effect = Observed effect")) +# scale_x_continuous(breaks = seq(from = 30, to = 100, by = 10), limits = c(30,100)) +labs(title ="Average Type M Error by sample size, effect size, and outcome standard deviation", y ="Type M error",x ="Sample Size",color ="Legend",caption ="Error bars represent 95% confidence intervals around the average type M error.") +theme_bw() +theme(plot.title =element_text(size =14),legend.title =element_text(size =14), legend.text =element_text(size =14),axis.title.x =element_text(size =14),axis.title.y =element_text(size =14),axis.text.x =element_text(size =12),axis.text.y =element_text(size =12),strip.text.x =element_text(size =14),plot.caption =element_text(size =12)) +facet_wrap(~ true_treat_eff + sd, ncol =2)mbias_plot
In the graph above, I plot the average ratio of the observed effect to the true effect (Type M error) conditional on observing a statistically significant result (p < .05). A ratio of 1 (black dotted line) indicates that the observed effect is equal on average to the true treatment effect.
Ratios bigger than 1 indicate that effects detected are exaggerated, ratios less than 1 indicate that the detected effects are smaller than the true effect, while a negative ratio indicates that the average detected effect has the opposite sign (direction) of the true simulated effect. To get a sense of the range of the ratio of observed to true effects, I included 95% confidence intervals around the average ratio.
As can be expected, lower sample sizes, smaller effect sizes, and a wider standard deviation of the outcome variable tend to increase the size and range of the Type M error and exaggerate the observed effect size. In the scenario with the lowest power (effect size of 3, outcome sd 10, and sample size of 30), the average Type M error is around 2.5 for the adjusted and 3.5 for the unadjusted models. The range of the Type M error is also incredibly large, spanning from -3 to 5 for the unadjusted model and -2 to 4.5 with adjustment, indicating that studies with very low power can produce statistically significant results in the opposite direction of the true effect - i.e. Type S or sign error (Gelman and Carlin 2014). Covariate adjustment also tends to push the size and range of the Type M error closer to the observed effect size (ratio goes to 1). As the treatment effect increases and the outcome standard deviation decreases, the Type M error for both the adjusted and unadjusted models tends to converge around 1, with a smaller range in the observed effect size.
References
Blair, Graeme, Alexander Coppock, and Macartan Humphreys. 2023. Research Design in the Social Sciences: Declaration, Diagnosis, and Redesign. Princeton. NJ: Princeton University Press. doi:10.7910/DVN/HYVPO5.
DellaVigna, Stefano, and Elizabeth Linos. 2022. “RCTs to Scale: Comprehensive Evidence From Two Nudge Units.”Econometrica 90(1): 81–116. doi:10.3982/ECTA18709.
Gelman, Andrew, and John Carlin. 2014. “Beyond Power Calculations: Assessing Type S (Sign) and Type M (Magnitude) Errors.”Perspectives on Psychological Science 9(6): 641–51. doi:10.1177/1745691614551642.
Simpson, Adrian. 2023. “A Recipe for Disappointment: Policy, Effect Size, and the Winner’s Curse.”Journal of Research on Educational Effectiveness 16(4): 643–62. doi:10.1080/19345747.2022.2066588.
